module HMGZ
!    real(8), pointer, dimension(:,:) :: phiftemp
    implicit none
contains

subroutine init_HMGZ
    include "pspec.FI"
    include "xsec.FI"
    
!    allocate (phiftemp(totmesh, ng))

end subroutine

subroutine sethk(k)
    include "hmgz.FI"
    real(8), intent(in)::k
    hk = k
end subroutine

subroutine homogenize_cycle
    use define
    use mvwnd
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: g, gp, i, i_f, i_c
    real(8) :: rtot, rrmv, rsct(ng), volflux, rvolflux
    
    ! flux condense
    hphi = 0.
    hphi2= 0.
    hTot=0;hRmv=0;hmSct=0;hNuF=0
    
!    ! cycle flux are used to obtain initial guess
!    do g=1, ng
!        do i=1, totmesh
!            i_c = fi2co(i)
!            hphi(i_c, g) = hphi(i_c, g) + phi_f(i,g)*h(i)
!        end do
!    end do
!    do g=1, ng
!        do i=1, totm_cmfd
!            hphi(i, g) = hphi(i, g)/h_cmfd(i)
!        end do
!    end do
    
    ! surface current condense   
    hpj = 0 
    do g=1, ng
        i_f = 1
        do i=1, totmesh*2
            if (jcvt2cmfd(i)) then
                hj(i_f, g) = jf(i, g)
                hpj(i_f, g) = pjf(i, g)
!                hphis(i_f,g) = phi_sa(i/2+1, g)
                i_f = i_f + 1
            end if
        end do
    end do

    do g=1, ng
!        rtot = 0;rrmv=0;rsct=0;
        do i=1, totmesh
            i_c = fi2co(i)

            volflux = phi_f(i,g)*h(i)
            hphi(i_c, g) = hphi(i_c, g) + volflux
            hTot(i_c,g) = hTot(i_c,g) + xTot(complist(i),g)*volflux
            hRmv(i_c,g) = hRmv(i_c,g) + xRmv(complist(i),g)*volflux
            hNuF(i_c,g) = hNuF(i_c,g) + xNuf(complist(i),g)*volflux
            do gp=1, ng
                hmSct(gp,g,i_c) = hmSct(gp,g,i_c) + XsctM(gp,g,complist(i))*volflux
            enddo 
            hphi2(i_c,g) = hphi2(i_c,g) + phi_f2(i,g)*h(i)/h_cmfd(i_c)
        end do
    end do
    do g=1, ng
        do i=1, totm_cmfd
            if (hphi(i, g)>0.) then
                rvolflux = 1./hphi(i, g)
                hTot(i,g)=hTot(i,g)*rvolflux
                hRmv(i,g)=hRmv(i,g)*rvolflux
                hNuF(i,g)=hNuF(i,g)*rvolflux
#ifdef _LANISO
            hDif(i,g)=1./(3*(hTot(i,g)-0.1_8))            
#else
            hDif(i,g)=1./(3*hTot(i,g))
#endif
                hDif(i,g)=1./(3*hTot(i,g))
                do gp=1, ng
                    hmSct(gp,g,i) = hmSct(gp,g,i)*rvolflux
                end do
                hphi(i, g) = hphi(i, g)/h_cmfd(i)
            endif
        end do
    end do
end subroutine

subroutine homogenize_acc
    use define
    use mvwnd
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: g, gp, i, i_f, i_c
    real(8) :: rtot, rrmv, rsct(ng), volflux, rvolflux
    real(8) :: vol(totm_cmfd)
    
    ! flux condense
    hphi = 0.
    hphi2 = 0.
    hTot=0;hRmv=0;hmSct=0;hNuF=0
    hTot_volavg=0;vol=0;
    gk = gkavg
    
!    ! cycle flux are used to obtain initial guess
!    do g=1, ng
!        do i=1, totmesh
!            i_c = fi2co(i)
!            hphi(i_c, g) = hphi(i_c, g) + phifAC(i,g)*h(i)
!        end do
!    end do
!    do g=1, ng
!        do i=1, totm_cmfd
!            hphi(i, g) = hphi(i, g)/h_cmfd(i)
!        end do
!    end do
    
    ! surface current condense    
    hpj = 0
    do g=1, ng
        i_f = 1
        do i=1, totmesh*2
            if (jcvt2cmfd(i)) then
                hj(i_f, g) = jfAC(i, g)
                hpj(i_f, g) = pjAC(i, g)
!                hphis(i_f,g) = phisaAC(i/2+1, g)
                hphis(i_f,g) = pjf(i/2+1, g)
                i_f = i_f + 1
            end if
        end do
    end do
    
    do g=1, ng
!        rtot = 0;rrmv=0;rsct=0;
        do i=1, totmesh
            i_c = fi2co(i)

            volflux = phifac(i,g)*h(i)
            hphi(i_c, g) = hphi(i_c, g) + volflux
            hTot(i_c,g) = hTot(i_c,g) + xTot(complist(i),g)*volflux
            hRmv(i_c,g) = hRmv(i_c,g) + xRmv(complist(i),g)*volflux
            hNuF(i_c,g) = hNuF(i_c,g) + xNuf(complist(i),g)*volflux
            do gp=1, ng
                hmSct(gp,g,i_c) = hmSct(gp,g,i_c) + XsctM(gp,g,complist(i))*volflux
            enddo 
            hphi2(i_c,g) = hphi2(i_c,g) + phi_f2(i,g)*h(i)/h_cmfd(i_c)
            
        end do
    end do
    do g=1, ng
        do i=1, totm_cmfd
            if (hphi(i, g)>0.) then
                rvolflux = 1./hphi(i, g)
                hTot(i,g)=hTot(i,g)*rvolflux
                hRmv(i,g)=hRmv(i,g)*rvolflux
                hNuF(i,g)=hNuF(i,g)*rvolflux
#ifdef _LANISO
            hDif(i,g)=1./(3*(hTot(i,g)-0.1_8))            
#else
            hDif(i,g)=1./(3*hTot(i,g))
#endif
                hDif(i,g)=1./(3*hTot(i,g))
                do gp=1, ng
                    hmSct(gp,g,i) = hmSct(gp,g,i)*rvolflux
                end do
                hphi(i, g) = hphi(i, g)/h_cmfd(i)
                
            endif
        end do
    end do
end subroutine

subroutine homogenize_acc_reduce(act, inact)
    use define
    use mvwnd
    implicit none
    include "pspec.FI"
    include "arrays.FI"
    include "accm.FI"
    include "xsec.FI"
    include "hmgz.FI"
    integer :: g, gp, i, i_f, i_c
    real(8) :: rtot, rrmv, rsct(ng), volflux, rvolflux
    integer, intent(in) :: act, inact
    
    ! flux condense
    hphi = 0.
    hphi2 = 0.
    hTot=0;hRmv=0;hmSct=0;hNuF=0
    gk = gkavg
    
!    ! cycle flux are used to obtain initial guess
!    do g=1, ng
!        do i=1, totmesh
!            i_c = fi2co(i)
!            hphi(i_c, g) = hphi(i_c, g) + phifAC(i,g)*h(i)
!        end do
!    end do
!    do g=1, ng
!        do i=1, totm_cmfd
!            hphi(i, g) = hphi(i, g)/h_cmfd(i)
!        end do
!    end do
    if (act .le. inact) then
        do g=1, ng
            phifac(:, g) = phifac(:,g) - phird(:,g,act)
            jfac(:, g) = jfac(:,g) - jrd(:,g,act)
            
        end do
    end if
    
    ! surface current condense    
    hpj = 0
    do g=1, ng
        i_f = 1
        do i=1, totmesh*2
            if (jcvt2cmfd(i)) then
                hj(i_f, g) = jfAC(i, g)
                hpj(i_f, g) = pjAC(i, g)
                hphis(i_f,g) = phisaAC(i/2+1, g)
                i_f = i_f + 1
            end if
        end do
    end do

    do g=1, ng
!        rtot = 0;rrmv=0;rsct=0;
        do i=1, totmesh
            i_c = fi2co(i)

            volflux = phifac(i,g)*h(i)
            hphi(i_c, g) = hphi(i_c, g) + volflux
            hTot(i_c,g) = hTot(i_c,g) + xTot(complist(i),g)*volflux
            hRmv(i_c,g) = hRmv(i_c,g) + xRmv(complist(i),g)*volflux
            hNuF(i_c,g) = hNuF(i_c,g) + xNuf(complist(i),g)*volflux
            do gp=1, ng
                hmSct(g,gp,i_c) = hmSct(g,gp,i_c) + XsctM(g,gp,complist(i))*volflux
            enddo            
            hphi2(i_c,g) = hphi2(i_c,g) + phif2AC(i,g)*h(i)/h_cmfd(i_c)
            hTot_volavg(i_c, g) = hTot(i_c,g) + xTot(complist(i),g)*h(i)
        end do
    end do
    do g=1, ng
        do i=1, totm_cmfd
            if (hphi(i, g)>0.) then
                rvolflux = 1./hphi(i, g)
                hTot(i,g)=hTot(i,g)*rvolflux
                hRmv(i,g)=hRmv(i,g)*rvolflux
                hNuF(i,g)=hNuF(i,g)*rvolflux
#ifdef _LANISO
            hDif(i,g)=1./(3*(hTot(i,g)-0.1_8))            
#else
            hDif(i,g)=1./(3*hTot(i,g))
#endif
                hDif(i,g)=1./(3*hTot(i,g))
                do gp=1, ng
                    hmSct(g,gp,i) = hmSct(g,gp,i)*rvolflux
                end do
                hphi(i, g) = hphi(i, g)/h_cmfd(i)
                hTot_volavg(i,g)=hTot_volavg(i,g)/h_cmfd(i)
            endif
        end do
    end do
end subroutine
subroutine homogenize_wnd

end subroutine
    
subroutine balance_control
    include "hmgz.FI"
    include "pspec.FI"
    include "xsec.FI"
    integer :: i, g, gp
    real(8) :: leak, rtot, rfis, rsct
    
    do g=1, ng
        do i=1, totm_cmfd
            leak=hj(2*i-1,g)+hj(2*i,g)
            rtot=htot(i,g)*hphi(i,g)
            rfis=0; rsct=0;
            do gp=1, ng
                rfis = rfis + hnuf(i,gp)*hphi(i,gp)
                rsct = rsct + hmsct(g,gp,i)*hphi(i,gp)
            end do
            hphi(i,g) = (rfis/hk+rsct-leak)/htot(i,g)
        end do
    end do
end subroutine    

    
end module
